Hard sphere crystallization gets rarer with increasing dimension 
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We recently found that crystallization of monodisperse hard spheres from the bulk fluid faces 
a much higher free energy barrier in four than in three dimensions at equivalent supersaturation, 
due to the increased geometrical frustration between the simplex-based fluid order and the crystal 
[J. A. van Meel, D. Frenkel, and P. Charbonneau, Phys. Rev. E 79 030201(R) (2009)]. Here, we 
analyze the microscopic contributions to the fluid-crystal interfacial free energy to understand how 
the barrier to crystallization changes with dimension. We find the barrier to grow with dimension 
and we identify the role of polydispersity in preventing crystal formation. The increased fluid 
stability allows us to study the jamming behavior in four, five, and six dimensions and compare our 
observations with two recent theories [C. Song, P. Wang, and H. A. Makse, Nature 453 629 (2008); 
G. Parisi and F. Zamponi, Rev. Mod. Phys, in press (2009)]. 

PACS numbers: 05.20.-y, 61.20.-p, 64.70.dm, 64.70.Q- 



Structural glasses form under conditions where, though 
the thermodynamically stable state of the system is crys- 
talline, the supersaturated fluid remains disordered. In- 
stead of nucleating crystals, the fluid phase becomes 
steadily more viscous, until the microscopic relaxation 
processes become slower than the experimental or simu- 
lation timescales. A glass is then obtained. To avoid en- 
countering a kinetic spinodal before falling out of equilib- 
rium, good glass formers should therefore be poor crys- 
tallizers [l|. Geometrical frustration is one of the fac- 
tors thought to slow the formation of ordered phases 
and thus assist glass formation |2|. Simple isotropic liq- 
uids are considered geometrically frustrated because the 
tetrahedron-based local order of the liquid cannot pack 
as a regular lattice. This scenario contrasts with what 
happens in a fluid of two-dimensional (2D) disks, where 
triangular order is locally as well as globally preferred 
and crystallization is particularly facile. 

The initial formulation of geometrical frustration by 
Frank considered the optimal way for kissing spheres in- 
teracting via a Lennard-Jones model potential to clus- 
ter around a central sphere Q. Frank found the icosa- 
hedral arrangement to be more energetically favorable 
than the cubic lattice unit cells. Though the original 
argument relies on the energetics of spurious surface ef- 
fects jH, mean- field solvation corrections leave the re- 
sult unchanged d @. For hard spheres, the icosahe- 
dron, which is the smallest maximum kissing-number 
Voronoi polyhedron, is also the optimally packed clus- 
ter. From the entropic standpoint the optimality of the 
structure therefore remains. More recently geometrical 
frustration has been couched in terms of the spatial cur- 
vature necessary to permit a defect-free lattice packing 
of tetrahedra (or, more generally, simplcxcs), which are 
the smallest building block in a Delaunay decomposition 
of space @, Q ■ This polytetrahedral scenario ascribes the 



presence of icosahedra to their singularly easy assembly 
from quasi-regular tetrahedra. Our recent study of 4D 
crystallization confirmed the simplex-based order in sim- 
ple fluids as the source of frustration. The 4D optimal 
kissing cluster, which can also tile space in the densest 
known lattice packing, plays however no identifiable role 
in the liquid order [9(. The observation that optimal 
clusters, such as the icosahedron, are not singular is in 
agreement with the careful examination of the local fluid 
structure [H El, an d offers a reasonable explanation 
for the limited amount of icosahedral order identified in 
experiments [U, El El El and simulations [M El of 
supercooled fluids. 

Geometrical frustration also contributes to the nucle- 
ation barrier. In the absence of impurities or inter- 
faces, crystallization proceeds through homogeneous nu- 
cleation in supersaturated solutions, as was spectacul arly 
observed in container- less levitated metallic liquids [18l j . 
According to classical nucleation theory (CNT), homo- 
geneous nucleation occurs through a rare fluctuation, 
whereby a crystallite that is sufficiently large for the bulk 
free energy gain to outweigh the interfacial free energy 
cost forms [19]. Crystallization then spontaneously pro- 
ceeds. The free energy difference between the ordered 
and disordered phases is fairly well understood micro- 
scopically in terms of the crystal packing efficiency. But 
the interfacial free energy contains a geometrical frus- 
tration contribution that is more challenging to inter- 
pret [11. Monodisperse hard spheres, the simplest sys- 
tem in which to study these effects, can indeed be super- 
cooled in 3D, but they are notoriously bad glass formers. 
Our earlier study, which found crystallization barriers 
in 4D to be much higher than in equivalently supersat- 
urated 3D fluids [9J, suggest that 3D hard spheres are 
only moderately geometrically frustrated. In this arti- 
cle we improve on this qualitative assessment: to quan- 



tify geometrical frustration, we look at the dimensional 
trends and use bond-order parameters and the fluid-hard- 
wall interfacial free energy as a reference. Equipped with 
a microscopic understanding of geometrical frustration, 
we examine some of its consequences for 3D polydis- 
perse spheres and consider how hard sphere crystalliza- 
tion evolves in higher dimensions. We find crystalliza- 
tion to become very rare, as was previously observed in 
systems of up to 6D [2l|. This rarity allows us to con- 
sider the consequences of the deeply supersaturated fluid 
branch on jamming, and compare the jamming results 
with the predictions of two recent theories [22|, [23| . 

The plan of the paper is as follows. First, we complete 
the phase diagram of 5D and 6D hard spheres, and use 
quantitative tools to describe the fluid and crystal orders 
in the various dimensions. We then compute the fluid- 
hard-wall interfacial free energy in 2D, 3D and 4D, in 
order to quantify the contribution of geometrical frustra- 
tion to the fluid-crystal interfacial free energy. Finally, 
we use these results to obtain the behavior of the free 
energy barrier to crystallization in higher dimensions. 



I. METHODOLOGY 

For convenience, in the rest of this article the particle 
diameter a sets the unit of length and the thermal energy 
ksT sets the unit of energy. For hard interactions this 
choice can be done without loss of generality, because 
entropy is the sole contributor to the free energy. 

Spheres become less efficient at filling space with in- 
creasing dimension. Though with our choice of units the 
fluid densities of interest p increase, the corresponding 
volume fraction 77 



77 = P V d 2- 



T d/2 



r(i + d/2)2 a 



(i) 



steadily decreases, because the volume of a d-dimensional 
hard sphere of radius 1/2 in K d , Vd2~ d , decreases faster 
than the crystal density increases. We report most quan- 
tities as volume fractions, but we also at times use p if it 
simplifies the notation. 



A. Phase Diagram 

Because of computational limitations, 6D is the maxi- 
mal dimension for which the phase diagram can reliably 
be obtained by simulation at this point. In a given dimen- 
sion, in addition to the fluid phase we consider the crystal 
phase postulated to be the densest and a less dense crys- 
tal for reference. The densest known close-packed struc- 
tures in 5D and 6D are degenerate through layering, the 
same way that hexagonal close-packed and face-centered 
cubic (fee) crystals are degenerate through layering in 
3D [24j |. For convenience we choose the most symmet- 
ric of these as order phases, which are D5 in 5D and 



Eq in 6D respectively [25| (see Appendix [5}. As in 3D, 
the impact of this decision on the phase diagram should 
be minimal [26j]. With increasing dimensionality layered 
structures show a growing similarity in their local two- 
and three-particle distributions, because layering affects 
only one of a growing number of spatial dimensions. The 
choice of specific layered phase should thus have but a 
small impact on the structural analysis. 

In order to precisely locate the freezing point, we com- 
pute the fluid and crystal hard sphere equations of state 
(EoS). Constant number of particles N, volume V, and 
temperature T Monte Carlo (MC) simulations [l?) give 
the radial pair distribution function g(r), which once ex- 
trapolated at contact is related to the EoS 



P/p= l + B 2 pg(l+), 



(2) 



where P is the pressure and B2 — Vd/2 is the second 
virial coefficient [28|. A sufficient number of MC cycles 
are used for the pressure to converge. Higher densities 
thus require longer simulations, because microscopic re- 
laxation becomes sluggish. A minimum of 50,000 MC 
cycles are performed, but up to ten times that amount 
is used when needed. Starting configurations are ob- 
tained by slowly compressing the system and by equi- 
librating at each density of interest along the way. For 
the fluid compressed beyond the freezing point, no crys- 
tallization is detected, which allows to thoroughly sample 
the metastable fluid, up until the microscopic relaxation 
becomes longer than the simulation time. The EoS are 
obtained for systems of 2048 (D 4 ), 4096 (4D fluids and 
A 4 ), 3888 (5D fluids and D 5 ), 14400 (A 5 ), 2048 (As), 
10000 (6D fluids), and 17496 (E 6 ) particles. To locate 
the fluid-crystal coexistence regime, we determine the 
absolute Helmholtz free energy per particle / of the crys- 
tal using the Einstein-crystal method [29] at a reference 
point: for D4 and A4 crystals we use 77 = 0.37, for As 
and 77 = 0.21, and for A and Eq 77 = 0.12. The excess 
free energy at other crystal densities is then obtained by 
thermodynamic integration of the EoS [13]. To obtain 
the fluid free energy the EoS is integrated from the ideal 
gas limit. The chemical potential 



p(p) = f + P/p 



(3) 



gives the fluid-crystal coexistence pressure p coex and 
^cocx ky finding where it is equal for both phases 
Ap(P) — 0. The densities of the coexisting phases is then 
obtained from constant N, p cocx ; and T MC simulations. 
This approach is formally equivalent to the common tan- 
gent construction, but we find it to be numerically more 
efficient. 



B. Order Parameters 

To characterize the structure of the fluid and crys- 
tal phases we need a criterion to quantify local order- 
ing. Studies in 2D and 3D suggest that order param- 
eters derived from rotationally-invariant combinations 
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of the m different spherical harmonics Y™ of degree I 
might suffice [13, HH, [Hj]. Here, we consider second- 
and third-order invariants, which are sensitive to the de- 
gree of spatial orientational correlation of the vectors 
that join neighboring particles. For a proper choice of 
I the (renormalized) invariant's (absolute) value is one 
for a perfect crystal and close to zero for a perfectly 
isotropic fluid. Though a 4D canonical spherical har- 
monics basis [H, Sec. 9.6] and both its second- and third- 
order invariants [34l ] are known, in higher dimensions it 
rapidly becomes analytically intractable to identify a ba- 
sis composed of weight vectors for the representation of 
SO(d) [35]. ft is therefore more convenient to rewrite 
the invariants as polynomials of the vector inner prod- 
ucts. For the second-order invariants, one simply uses 



the Gegenbauer polynomials G 



d/2-l 

l 



obtained from the 



sum rule [36l . Thm. 9.6.3]. For instance, the sum over 
the (I + l) 2 4D spherical harmonics for unit vectors 
can be rewritten as 



2tt 2 



O+i) 2 

]r Yr(h)Yr(h). (4) 



m— 1 



The second-order local bond-order correlator qi(i,j) is 
then obtained by summing over the N(i) and N(j) near- 
est neighbors of particles i and j conveniently chosen to 
be within a distance equal to the first minimum of g(r). 
By letting the indices a and (3 run over these neighbors 
we find [9j 



v JV(i) v W(i) r d/2-i r „ . , 
N{i)N(j) 

This two-particle second-order invariant qi(i,j), rather 
than a single-particle qi(i,i) or a global second-order in- 

variant Qi = N(i)N(j)qt(i, j)] * * /£. N(i), is 



known to be a reliable basis for the identification of in- 
dividual crystalline particles in 3D [3l| and 4D [9|], and 
allows for a more sensitive analysis of the fluid structure. 

We now develop an approach to obtain third-order 
rotationally invariant polynomials wi analogous to the 
Gegenbauer polynomials. This approach, like the one 
described above for the second-order invariant, has the 
important advantage that we do not need a prior knowl- 
edge of a canonical spherical harmonics basis. A classical 
theorem due to Weyl says that any polynomial in m sets 



of variables X\, 
nal action 



. , X m 6 



invariant under the diago- 



g ■ f(Xt, ...,X m ) = f(gX u gX m ) for g E SO(d) 

can be expressed in terms of the inner products (X^, Xj) 
and the determinants detfX^ ...X id ]. For third-order 
invariants (m = 3) in d > 4, all the determinants are zero, 
and we are able to write the invariant polynomial in X — 
(X u ...,X d ), Y = (Yi, . . . , Yd), and Z - (Z 1 ,...,Z d ) 
in terms of the various inner products x = (X, X) , y — 
(Y,Y),z = (Z, Z) anda= (X,Y),b= (X,Z),c = (Y,Z). 

Let / be a polynomial in X, Y, Z. Suppose that / is 
invariant under the diagonal action of SO(e?). Then there 
is a polynomial p(x, y, z, a, b, c) such that 

f(X, Y, Z) = p(x, y, z, a, b, c). 

Lemma 1.1 Suppose f is homogeneous of degree I in X, 
Y , and Z separately, and is therefore homogeneous of 
degree 31 overall. Then 

p(X 2 x, fj?y, v 2 z, A/ia, Xvb, five) — {\pLi>) l p(x, y, z, a, 6, c). 
In particular, p is homogeneous of degree 31/2 overall. 
Let 



D x {p) = 2d 
D Y (p) = 2d 
D z {p) 



dp 

dx 
dp 
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4c 



4c 
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d 2 p 
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d 2 p 
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d 2 p 
. — — 4 
dbdc 
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d 2 p 
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db 2 ' 
d 2 p 
dc 2 ' 
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9c 2 ' 



Lemma 1.2 The operators above satisfy 

^\p{x,y,z,a,b,c) = (D x p)(x,y,z,a,b,c), 

and similarly for Y and Z. In particular if f is harmonic 
in X, Y and Z, then D x {p) = D Y (p) = D z (p)=0. 

The proof is an exercise in using the chain rule. 



Using Lemmas 11.11 and 11.21 we set up a system 
of equations for the coefficients of the polynomial 
p(x, y, z, a, 6, c) corresponding to a SO(d) invariant poly- 
nomial f(X, Y, Z) of degree I and harmonic in each of 
X, Y, and Z separately. We can solve this system of 
equation, and once we choose the normalization, say 



3 



p(l, 1,1,0, 0,0) = 1, there is a unique solution. Set- 
ting x = y — z = 1 restricts the obtained func- 
tion on the unit sphere. We call the resulting func- 
tion wf(a,b, c). Examples are given in Appendix IBl For 
the reader cognisant of representation theory, note that 
wf(X,Y,Z) = wf(a,b,c) is a basis for the copy of the 
irreducible representation in the triple tensor product 
HiiS*- 1 ) ® HiiS^ 1 ) <g> HiiS^ 1 ). 

As in the second-order case, those polynomials allow 
us to rewrite the third-order local bond-order correlator 
Wi(i) up to a dimension-dependent multiplicative con- 
stant c. 



the coupling parameter A = oo, hard spheres are con- 
fined by a hard wall, while for A = the bulk fluid limit 
is recovered. The interfacial free energy is obtained by 
Kirkwood integration 



Tf- 



1 

2A 



dX 



dX 



(9) 



where A the area of a single side of the wall. In practice 
the integral is solved using a 21 point Gaussian-Kronrod 
formula in a finite interval A £ (0, A max ) with A max chosen 
to approximate arbitrarily closely a hard wall. 



a y s ' r »/3. r »q ■ TiS,r il3 ■ v lS ) 

T 2"-2[JV(i)]3[q,(i,i)]3/2 ' 



Wi(i) 



(6) 



In 3D and 4D, the constant cf can be set by compar- 
ing with the expression available in the literature. In 
higher dimensions, we choose the normalization for which 
the polynomial equals unity when evaluated on three or- 
thogonal unit vectors i.e., cfwf(0, 0,0) = 1. Note that 
because of the rotational symmetry, triplets with per- 
muted indices can be summed only once by correcting 
for the multiplicity. This simplification offers an im- 
portant computational advantage. Though the use of 
rotationally-invariant polynomials for the computation 
of the bond-order parameters is mainly used for analyt- 
ical convenience, it is also worth noting that for large 
I and at low densities, it is computationally more effi- 
cient than the standard spherical harmonics decomposi- 
tion, and that their algebraic simplicity minimizes the 
risks of error at the implementation stage. 



C. Wall Cleaving Surface Tension 

The 2D, 3D, and 4D hard sphere fluid-hard-wall inter- 
facial free energy 7f_ w is calculated through MC simu- 
lation using the higher-dimensional generalization of an 
earlier thermodynamic integration scheme [37j . We start 
from a system periodically confined by both sides of a 
hard wall perpendicular to the x axis and make the wall 
gradually more penetrable until the bulk system is ob- 
tained. Confinement is achieved by introducing the aux- 
iliary Hamiltonian 



A' 



N 



(7) 



hi 



i = l 



where Vhs is the hard-core exclusion between hard 
spheres and V w is the penetrable wall potential 



V w (x) 



exp(-2x) if \x\ < er/2 
otherwise 



(8) 



for a sphere a distance x from the wall. This truncated 
exponential function is known to provide a high numeri- 
cal accuracy route for 7f _ w computation [37l [3q | . When 



D. Generalized Classical Nucleation Theory 

Classical nucleation theory (CNT) [l9| considers con- 
tributions from chemical potential difference between the 
bulk phases and the fluid-crystal interfacial free energy 
7f _ x of a spherical crystallite to obtain a free energy func- 
tional 

AG(n) = A e2 (n/p x ) (d - 1)/d 7f-x - nAfi, (10) 

of the number of particles n in the crystallite. The func- 
tional further depends on the crystal density p x at the 
supersaturated fluid pressure and on a geometrical pref- 
actor Ad- For hard spheres Ad = Sd-iV^ d 1 , where 
Sd—i — dVd is the surface area of a rf-dimensional unit 
sphere. The resulting barrier height at the critical cluster 
size n* is 



AG*(n*) = 



it, 



r(d/2 + l)2<* pi-i 



(11) 



and in the high dimensional limit the barrier asymptoti- 
cally approaches 



AG* (71*) - (27red) 



d/2_ 



7f- 



1 An 



d-l 



(12) 



The rate of nucleation per unit volume k is then 
k = Kcxp(— AG*), where k is a kinetic prefactor propor- 
tional to the diffusion coefficient in the fluid phase 32l |. 
The kinetic prefactor has a weak dimensionality depen- 
dence that we won't consider here. Though schematic, 
this level of theory is sufficient to clarify the contribution 
of geometrical frustration to the crystallization barrier 
through an analysis of 7f_ x . Within the CNT framework 
geometrical frustration between ordered and disordered 
phases should lead to a relatively large 7f _ x , and thus to a 
high crystallization free energy barrier. On the contrary, 
geometrically similar phases should have small 7f_ x and 
AG*(n*). 



II. RESULTS AND DISCUSSION 
A. Phase Diagram and Jamming 

The computed fluid EoS agrees with earlier 4D [§], 
5D [H, EHH, Ei|, and 6D g| simulation results as 
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FIG. 1: (Color online) Monte Carlo EoS of the fluid and the two densest known ordered phases in 4D [9[ (a), 5D (c), and 6D 
(d) computed at constant V (filled symbols) and constant P (open symbols), with Pade approximants of the virial expansion 
for the fluid [2||, and the Speedy fits to the crystal phase results [3^ (solid lines). Insets give A/i and the common tangent 
construction for determining coexistence between the fluid and the densest crystal phase. The additional panel for 4D (b) 
contrasts MC and SPT/CT EoS as well as the resulting coexistence determination (inset). 
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TABLE I: Coexistence parameters from Monte Carlo simulations compared with previous simulation estimates (see text) and 
the corrected SPT/CT results (see text). The volume fraction of the densest known lattice ry cp is also included for reference [25l . 
Chap 1. § 1.5] 



well as the 5-4 Pade approximants of the virial expan- 
sion [H, EH (Fig. Q]) . Small deviations are only observed 
at the highest densities, where the expansion is less ac- 
curate [28l [45| . Crystal phase EoS for the D4 and -D5 
lattice geometries in 4D and 5D respectively were first 
obtained from simulation in the early 80 's, but without 
reference free energies [U |42j , and we are not aware of 
any 6D simulation results. As expected from free volume 
arguments, the densest known lattice is the phase with 



the lowest pressure at densities where it is mechanically 
stable, and is the most free energetically favorable of the 
ordered phases. Assuming that the crystallization ki- 
netics is controlled by the free energy barrier height, the 
most stable ordered phase should be the only relevant one 
of hard spheres. The generation of crystallites with other 
symmetries is only possible at much higher pressures and 
with a smaller thermodynamic drive. The scaling of the 
fluid-crystal interfacial free energy with pressure energy 
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TABLE II: Volume fraction of the maximally-random jammed 
states obtained from free- volume extrapolation 7/mrj, direct 
compression [2l| rfj°^ , an extension of the statistical mean- 
field approach of Ref. [22|] ??mrj j an d the replica method mean- 
field approach [23| 



suggests that neither phase would have a significant ad- 
vantage over the other from that respect (cf. Sect. Ill Dj) . 

The fluid-crystal coexistence conditions from simula- 
tion for hard spheres in 5D and 6D, along with the ear- 
lier 3D [4l| and 4D [9J] results are reported in Table fl] 
Skoge et al. offered upper bounds to the 4D and 5D co- 
existence regimes by using the pressure of the fluid at 
the density at which the simulated crystal becomes 
mechanically unstable as an estimate of coexistence pres- 
sure p coex [2l|. A more accurate estimate of p cocx can 
be obtained from the same data by using instead a quasi- 
Maxwell construction around the limit of mechanical sta- 
bility [47( . We include the results of this last analysis and 
the coupled fluid scaled-particle theory (SPT) and crys- 
tal cell theory (CT) coexistence determination of Finken 
et al. [H in Table[flas well. To the best of our knowledge, 
density functional theory (DFT) coexistence parameters 
have only been reported for the non-equilibrium 4D fluid- 
A4 pair of phases [48[, which does not lend itself to a 
meaningful comparison with simulations. Finken et al. 
do refer to DFT coexistence calculations, but do not re- 
port their results for 4D to 6D do[. 

Our MC results are at least an order or magnitude 
more precise than the estimates from Ref. [2l|, which 
permits a clearer assessment of the SPT/CT predictions. 
The SPT/CT analysis correctly captures certain dimen- 
sional trends. For instance, the relative difference in crys- 
tal and fluid volume fraction at coexistence Ai] coox /rh, 
which is thought to go to unity for large dimensions [401 ] , 
does increase appreciably from below 10% in 3D to over 
20% in 6D. And the crystal volume fraction at coexis- 
tence rj x decreases relatively faster than the close-packed 
volume fraction rj cp , which leaves the phase diagram dom- 
inated, in percentage of the accessible densities, by the 
ordered phase [49]. Finken et al., however, incorrectly 
implement the common tangent construction, which ex- 
plains why their reported coexistence regimes [4oT ] are 
well above our simulation results. Once corrected, the 
theoretical SPT/CT predictions more closely follow the 
simulation results (see Table|T]). SPT shows a fairly good 
agreement with the fluid EoS, and, as expected from the 
third-order virial coefficient, overshoots the fluid pressure 
at high densities [4(| (see Fig. [T]). CT however signifi- 
cantly underestimates the crystal pressure near coexis- 
tence and the effect does not go away with dimension. 
The high compressibility of hard sphere crystals near the 



limit of mechanical stability is a collective effect that is 
not captured by the mean-field nature of the theory. The 
cancelation of errors leaves SPT/CT coexistence densi- 
ties reasonably on target, but SPT/CT overestimates the 
width of the coexistence densities and the coexistence 
pressure. It is unclear if this effect vanishes with dimen- 
sion. 

It is interesting to note that p cocx does not change 
monotonically with dimension, but goes through a mini- 
mum in 4D. The nonmonotonic behavior of p cocx might 
be due to D^s particularly well-suited nature to fill 4D 
Euclidian space. A D4 lattice can be generated by placing 
a sphere in each of the voids of a 4D simple cubic lattice. 
These new spheres are equidistant to the ones on the sim- 
ple cubic frame, so the resulting lattice is twice denser. 
The corresponding construction P3, or body centered- 
cubic lattice, in 3D packs much less efficiently, because 
the simple cubic frame needs to be extended to insert 
the new spheres. Though D4 does not appear singular 
in the dimensional trend of dense packings [25|, Chap. 
1, §1.5], its specificity might simply be overshadowed by 
other dimensional trends to which p cocx is less sensitive. 
We therefore conjecture that the non-monotonic coexis- 
tence pressure is a symmetry signature that should also 
be present in 8D, 12D, 16D, and 24D, where other sin- 
gularly dense lattices are known to exist. 

Another noteworthy observation concerns the high 
pressure limit of the fluid EoS. The fluid results pre- 
sented in this article are only for systems in equilibrium 
or in metastable equilibrium, i.e. the initial configura- 
tions are prepared in the fluid state and the simulations 
are run much longer than the fluid microscopic relax- 
ation timescale, though much shorter than the nucleation 
timescale. But because no crystallization takes place on 
the simulation time scale, only the growing microscopic 
relaxation timescale limits the range of densities can be 
reached in simulations. The EoS can thus obtained for 
supersaturated systems at much higher pressures than in 
3D (see Sect. Ill E|) . We extract the infinite-pressure com- 
pression limit of the densest supersaturated fluid point on 
the EoS from the free- volume functional form for pressure 
suggested in Ref. [H 



V d 2* l-{ r] /^) 1,d 



(13) 



This functional form, which is supported by the analy- 
sis of Parisi and Zamponi [23| for 4D fluids, corresponds 
to a non-equilibrium compression so rapid that no mi- 
croscopic relaxation can take place. Because microscopic 
relaxation becomes increasingly slow large at high fluid 
densities, it approximates the compression algorithm of 
Ref. [2l|]. The volume fraction of the disordered jammed 
system corresponding to the densest equilibrated fluid is 
obtained by solving for rf^ . The parameter if^ is not the 
volume fraction of the maximally-dense random jammed 
(MRJ) state per se. But were it possible for the fluid to 
avoid crystallization indefinitely, which is obviously an 
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tial correlations vanish with increasing system dimen- 
sion [2111 . In accordance to the "decorrelation princi- 
ple" [51(, we also expect orientational correlations of or- 
der I, gi(r), to decay more rapidly when the dimension in- 
creases. As seen in Figure [2j in 2D the hexatic signature 
gives rise to long-ranged orientational correlations with a 
length that diverges on approaching coexistence [12, HH ; 
in 3D the orientational order stretches over a couple of 
particle radii, but stays finite even in the supersaturated 
regime; and in higher dimensions the correlations keep 
decaying, as can be assessed from the height ratio of the 
second to the first peak of ge (r) in Fig. O Note that a 
similar behavior is observed for all I considered, but I = 6 
has the advantage of capturing the crystal symmetry for 
all dimensions. 



FIG. 2: (Color online) Impact of dimensionality on the decay 
of orientational order. Top: Second to first peak ratio of 
the radial decay of the orientational order parameter ge(r) 
at fluid coexistence. The line is a guide to the eye. Bottom: 
Decay of ge(r) in the fluid near the hexatic phase in 2D and at 
coexistence in 3D and 4D. The 5D and 6D plots (not shown) 
are qualitatively similar to the 4D results. 



unphysica! but useful abstraction, it should asymptoti- 
cally approach it from below. A quenched compression 
started from the non-crystallizing, but otherwise equili- 
brated fluid branch at higher supersaturations, where mi- 
croscopic relaxation is even slower, would give a denser 
rf^ [23|. The high densities achieved here should thus 
give a fairly reliable lower bound to ?7mrj ■ The resulting 
values are in very good agreement with the volume frac- 
tions obtained by direct non-equilibrium compression [21( 
(see Table [TTJ) . A recent statistical mean- field theory of 
jamming that shows a surprisingly high accuracy with 
the experimental and simulation j/mrj in 3D [22j |. would 
be expected to perform similarly well, if not better, when 
dimensionality is increased, due to the growing number 
of nearest neighbors. On reproducing the arguments of 
Ref. [22j for higher dimensions (see Appendix [U]), we find 
that though the accuracy is still quite good, it is not 
as high as for 3D. The propagated relative error of the 
statistical mean-field treatment is about 5-10%, which 
suggests that the high accuracy of the analytical 3D re- 
sults of Ref. [22j might be fortuitous or that it is par- 
ticularly sensitive to the choice of scaling assumptions. 
The mean-field treatment based on the replica method 
offers however predictions that more consistent with the 
simulation results [23[ . The dimensional trends are more 
similar, and the simulation results are slightly smaller 
than the theoretical predictions, which is precisely where 
they are expected to be p3| . 



The authors of Ref. [2l| remarked that the number 
of particles counted in the first peak of g(r) for super- 
saturated fluids matches the number of kissing neigh- 
bors in the densest known lattice phase for a given di- 
mension. They hypothesized that disordered packings in 
higher dimension might thus be built of deformed crystal 
unit cells, in contrast to the three-dimensional case where 
"icosahedral" order was thought to dominate the pack- 
ing. The distributions of local bond-order correlators, 
which shows how the relative crystal and fluid local or- 
ders evolve with dimensionality paint a different picture 
(see Fig. Both the second- and third-order invari- 
ants in 4D to 6D, capture no significant overlap between 
the liquid and the crystal local bond-order parameters, 
in contrast to 2D and 3D where the bond-order distribu- 



B. Bond-Order Correlators 

Skoge et ah, considering the radial pair distribution 
function g(r), found that higher-order unconstrained spa- 



tions overlap significantly 5J|. Moreover, the distinction 
between the fluid and crystal phases increases with di- 
mension, which suggest that the fluid and crystal local 
orders are just as or more distinct with increasing di- 
mensionality, not less [5511 . In 4D, where the maximally- 
kissing cluster 24-cell [56| is also the unit cell of the crys- 
tal, but is not a simplex-based structure, no hint of the 
crystal order in the fluid is captured by the bond-order 
distribution (Fig. |3J) 0. For a simplex-based cluster to 
have as many nearest neighbors within the first peak of 
g(r) as in the crystal, the first neighbor spheres cannot 
all be kissing the central sphere at the same time, but 
have to fluctuate in and out of the surface of the central 
sphere. This variety of possible configurations is what 
broadens the firstpeak of g(r) and by ricochet its sec- 
ond peak as well [2lJ. Though they are harder to illus- 
trate geometrically, similar phenomena are expected in 
higher dimensions. The bond-order distribution is there- 
fore fully consistent with a fluid structure dominated by 
simplex-based order, but not with the presence of de- 
formed crystal unit cells. This clear separation in local 
order also suggests that for dimensions greater than three 
fluid configurations should be easier to distinguish from 
the partially crystalline or polycrystalline systems that 
can be observed in 3D compression studies [57] . 
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FIG. 4: (Color online) Top: Additional fluid volume per unit 
area of surface Ava = Jt- W /P from Monte Carlo simulations 
and the Pade approximant equation of state [28l . |58| ] (points) 
compared with the virial expansion of Eg, 1161 (lines). Bottom: 
Comparison of 7 X _ W /P from published simulation data [59l . 
IdQ| | and the Speedy equation of state [39( ] (points) with the 
cell-like theory (see text) for different orientations of the 3D 
fee crystals (lines). 



7f-x- We will get back to this point in Section III D I We 
first consider the behavior of the fluid-hard-wall interfa- 
cial free energy 7f_ w , which is easier to interpret micro- 
scopically and is a limit case for the fluid-crystal geomet- 
rical frustration in all dimensions. From the fluid point 
of view, both the hard wall and the crystal plane exclude 
configurations containing spheres less than a radius away 
from the hard surface, but the crystal has additional free 
volume crevasses to explore. 

The fluid-hard-wall interfacial free energy is an equi- 
librium quantity that, unlike the liquid-crystal interfacial 
free energy, is well defined at all densities where the liq- 
uid is stable. An ideal gas has no fluid-hard-wall inter- 
facial free energy, but the low density limit of spherical 
particles that exclude volume is 7f_ w = P/2 + 0{ P 2 ) = 
p/2 + 0(p 2 ), because of the PV work required to exclude 
particles from the surface of the hard wall. At higher den- 
sities, just like a surfactant reduces the surface tension 
by occupying part of the interfacial area, the presence of 
particles at the interface reduce the entropy cost for the 
other particles and partly offsets the increase in opposing 
bulk pressure. An exact expansion gives [6l|, [62|, [63| 



7f- 



P{P) 



- B mP 2 - B m p 3 - B mP 4 - 0{p 5 ), (14) 



C. Fluid-Hard- Wall Interfacial Free Energy 

The local bond-order distributions indicate that the 
hard-sphere fluid order resembles more the crystal order 
in 3D than in higher dimensions. Because of the purely 
entropic nature of hard-sphere systems, microscopic ge- 
ometry should have a clear thermodynamic signature on 
quantities such as the fluid-crystal interfacial free energy 



where the first couple of surface virial coefficients B$u 
are computed for hard spheres in all dimensions in Ap- 
pendix [D] 

The virial expansion can be contrasted with the SPT 
and the "mechanical" expressions for 7f_ w . It has al- 
ready been noted that 3D SPT, captures the first term 
in the virial expansion exactly and is within a few percent 
of the next order term [63|. This seems to be the case 
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for all odd dimensions. For even dimensions, though SPT 
captures the low-density pressure behavior correctly, it is 
off already for B n2 - In 2D SPT gives B n2 = n/8, while 
the correct value is 1/3. Even when using a precise ex- 
pression for the pressure, SPT docs not capture the inver- 
sion of 7 f c °°^ from 3D to 4D (see Table Ell). Similarly, it 
has also been pointed out that the "mechanical" result of 
Kirkwood and Buff with a Fowler-type approximation for 
the two-point correlation l65j (see also Ref. [H, Sect. 
4.4]) correctly captures the first surface virial coefficient 
in 3D [67| . and in all dimensions generally. Higher order 
terms contain surface effects that the simple decorrelated 
approximation cannot capture, which explains why it is 
about 60% off for Bq 3 in 3D, predicting _Bq 3 = 57r 2 /96 
though the exact value is 1497r 2 /1680. One should there- 
fore take quantitative 7f_ w predictions of those two ap- 
proximations with a grain of salt. 

To better understand the density dependence of 7f_ w , 
we remove its trivial pressure contribution and obtain the 
change in fluid volume per unit of wall area 



Av A (p) 



7f-w(p) 
P(P) ■ 



(15) 



Using the standard pressure virial coefficients Bi , we ex- 
pand Ava to second order in density 

Av A {p) = -- B m p + (B 2 Bn2 - B m )p 2 + 0{p 3 ). (16) 

For 3D hard spheres -Bn4 has been evaluated numeri- 
cally [63], so the expansion can be carried to the next 
order. As expected, the initial slope of Ava(p) is al- 
ways negative (see Table fVT) . Also, the coefficient of p 2 , 
which is initially negative, changes sign between 3D and 
4D. Because the p 2 coefficient is small in 3D and 4D, 
the next order coefficient is more significant and indeed 
Bq4 markedly improves agreement with 3D MC results. 
Fig. 2] shows that theory and simulations results nearly 
coincide at low density, which validates the current sim- 
ulation methodology and the latest 3D results [H^, [Hj] . 
For 3D the match with theory extends close to the co- 
existence limit. Only one or two additional coefficients 
would probably be needed to capture the full curve with 
simulation accuracy at coexistence. 

Based on the physical interpretation of Ava, we ex- 
pect the change in fluid volume to monotonically decay 
to a positive constant at high densities. By breaking the 
translational symmetry the introduction of a hard wall 
must perturb the fluid order to some extent. We also 
expect a plateau to develop on approaching that limit, 
because the high-density structure of the supersaturated 
fluid varies little. But it might only be possible to ob- 
serve this plateau at high fluid density, in the supersat- 
urated regime. In 2D high density reliable results are 
technically difficult to obtain. The possible presence of 
an hexatic phase at high densities commands very large 
system sizes, because the wall favors the local crystal or- 
der and pushes defects away from the interface. In 3D 




FIG. 5: Additional volume required for placing a wall through 
a simplex. The disk on the left must be pushed back suffi- 
ciently (arrow) to allow for the wall insertion (vertical line). 



rapid wall-induced crystallization occurs before a signifi- 
cant degree of saturation a reached [37|, |59j, [6CJ, [68j, [69( . In 
d > 3 however the hard interface does not accommodate 
a regular packing of simplexes, unlike in lower dimen- 
sions. We therefore expect heterogeneous crystallization 
to be sufficiently slow to study 7f_ w in the dense fluid 
regime by simulation. Though 4D Ava has not yet satu- 
rated in Fig. [H the results do suggest a slowdown of the 
decrease. Unfortunately, the higher-dimensional system 
sizes required to further clarify this point are currently 
beyond our computational reach. 

Indirect evidence nonetheless supports the saturation 
scenario. For crystals, saturation of Ava is similar to 
the cell theory that was successfully applied in 3D for 
7f-w [60( ] . The gap that must be opened to insert a plane 
along a given cut through a packed crystal is the high 
density limit of 2 Ava, where the factor of two accounts 
for the two interfaces that are created by wall insertion. 
Published crystal- wall interfacial free energy 7 X _ W simu- 
lation results for various faces of the fee crystal (5f| Irioj 
allow to extract the corresponding Ava{p)- The lower 
panel of Fig. 2] shows that, within the error bars, the 
simulation results follow pretty closely this simple satu- 
ration scaling. The toy model also rationalizes the "bro- 
ken bond" and anisotropy ratios of Ref. [6(| ■ Because the 
interfacial free energy between the fee (111) plane and a 
hard wall is tiled exclusively with 2D simplexes, which 
is similar to a 3D high density fluid-hard-wall interface, 
a similar Ava saturation is expected. This prediction 
can be checked by reanalyzing the recent simulations of 
jamming in confinement [70 ]. Desmond and Weeks con- 
sidered how the close-packed density of bidisperse spheres 
is affected when the distance between two confining hard 
walls is changed. The scaling quantity C they extract can 
be reexpressed as the additional fluid volume per wall sur- 
face area Ava(pmrj) — C?7mrj/2.8, where the numerical 
prefactor takes into account the presence of two interfaces 
and rescales the diameter of the larger spheres to unity. 
In 3D this gives Ava{pmkj) ~ 0.054, which is lower than 
the fee (111) limit for monodisperse spheres of 0.091. It 



9 



(1 


7f°x x 


cocx 
If — w 


Lj - L ^ 1 It — w 


SPTo-y? 00 * 

Ji 1 ^ It— w 


-Yf COCX /-V, COCX 
If— x / It— w 


3 


0.557 [59] 


1.98(1) 


2.30 [60] 


1.75 


0.28 


4 


1.0 


1.96(2) 


2.53 [40] 


1.89 


0.53 


5 






2.93 


2.07 




6 






4.05 


3.24 





TABLE III: Fluid-crystal and fluid-hard- wall interfacial free 
energies at coexistence from simulations and SPT [40l . [60| . 
The column SPTi is the SPT interfacial tension with SPT 
pressure, and the column SPT2 is the SPT interfacial tension 
calculated using the MC pressure. 



is qualitatively expected that a bidisperse system have a 
reduced Ava, because the smaller spheres do not need 
to be pushed as far back from the interface as the larger 
spheres [lllll. 

How the high-density limit of Ava changes with dimen- 
sion dictates whether the free energy barrier to nucleation 
vanishes or not in high dimensions. From the geometri- 
cal frustration analysis, we have identified the simplex as 
the dominant geometrical structure in high density fluids. 
It is therefore reasonable to expect that the dominant 
structure near a hard interface be a truncated simplex. 
Obviously the fluid interface does not exclusively contain 
truncated simplexes, so we further assume that simplexes 
are representative of the interface, and thus dominate 
the volume loss. The high density limit of Ava is then 
the distance that a vertex sphere needs to be pushed 
out from the other spheres in the simplex for a tangent 
plane to be inserted. For 2D, this scenario is depicted 
in Fig. \5\ and the details of the general calculation are 
given in Appendix [EJ Within this approximation Ava 
monotonically increases with dimension, but is bounded 
by 1/2 — \/2/4 ~ 0.146. In high dimensions, 7™^ does 
not therefore asymptotically vanish, but increases, due 
to the growing p coex . 



D. Fluid-Crystal Interfacial Free Energy 

Our previous crystallization study of 4D hard spheres 
gives values for 7f_ x that increase with fluid supersatu- 
ration (or, equivalently, pressure) @. These interfacial 
free energies are more than twice as large as for 3D hard 
spheres at comparable supersaturation 



32j. Even tak- 



ing into account the slightly different interfacial densities, 
the gap is large, justifying our earlier claim of increased 
geometrical frustration in 4D @. To allow for a more 
quantitative comparison, we consider here the system at 
coexistence, where the equilibrium fluid-crystal interfa- 
cial free energy 7™™ is unambiguously defined. 

Crystallization-derived 7f_ x for 3D hard spheres have 
been corrected for finite-size effects to obtain 7™^ x [z3- 
The effective critical cluster sizes contain at most a few 
spherical shells and the resulting strongly curved crys- 
tallite leads to a large internal Laplace pressure. The re- 
sulting loss of free volume per particle is compounded by 



the relatively large compressibility of higher-dimensional 
crystals near coexistence, increasing the free-energy cost 
of forming the interface. Just like 7f- w increases with 
density, however, the non-equilibrium fluid-crystal in- 
terfacial free energy for supersaturated fluids should be 
larger than jf °^ x , due to the overall increase in bulk pres- 
sure. Because Ava has probably not yet saturated, it is 
difficult to get a precise estimate of this effect, but it 
is of the right magnitude to explain the change of 7f_ x 
with pressure in both 3D and 4D. Which of the crystal- 
lite finite-size or the increase in fluid pressure dominates 
the non-equilibrium 'yf _ x behavior cannot be resolved 
here. But to first order, through the use of the Tolman 
"ansatz", they both suggest that the interfacial free en- 
ergy depends linearly on supersaturation A/i 74]. This 
scaling, though microscopically inaccurate and therefore 
only used on an ad hoc basis [75[ , showed a relative suc- 
cess in 3D [Zl, [ll, and gives 7 f c ° e x x w 1.0 for 4D hard 
spheres. In spite of its crudeness, the result is sufficiently 
precise to interpret the thermodynamic consequences of 
geometrical frustration, because the equivalent quantities 
in 3D are known with high accuracy 



quivalen 
59, 77]. 



To a first approximation one expects the fluid-crystal 
interfacial free energy to scale linearly with the fluid- 
hard-wall quantity if geometrical frustration were con- 
stant in all dimensions, because the depth of the inter- 
face remains of the order of the particle dimension. Yet 
comparing the ratio of fluid-crystal to fluid-hard-wall in- 
terfacial free energies in 3D and 4D in Table IIIII shows 
this not to be universally the case. To understand the ori- 
gin of the difference we consider how the interfacial free 
energy decreases upon changing the hard wall to a crys- 
tal interface. The core of the reduction comes from two 
sources: the crystal planes near the fluid interface have 
more free volume than those in the bulk and the interfac- 
ing fluid requires less Ava than next to hard wall. The 
interfacial picture suggested by DFT, where in 3D the 
bulk of the interfacial free energy comes from the set of 
three or four layers that form the interface, is consistent 
with both scenarios [78(- In 3D, the crystal relaxation 
due to the additional free volume for interfacial particles 
has recently been considered the dominant effect [79| . 
But an older toy model by Spaepen ascribes a signifi- 
cant contribution to the fluid [80l[8ll|. Though our results 
cannot resolve quantitatively the relative contribution of 
the two scenarios, they at least indicate that both are 
of comparable magnitude [13] • If the crystal modulation 
alone explained the reduction of the interfacial free en- 
ergy from the hard wall limit, it should be similar if not 
greater in 4D than in 3D, because the reduced p cocx in 
4D allows for more wandering of the interfacial crystalline 
particles. Yet the fluid-crystal interfacial free energy is 
both proportionally and absolutely larger in 4D than in 
3D. The further reduction of 7f ° x x in 3D must therefore 
arise from the relatively good geometrical match between 
the fluid and the crystal order at the interface. The crys- 
tal surface allows for more additional microstates in the 
spacing between the interfacial crystal spheres in the 3D 
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fee lattice than in the 4D D4 lattice. Based on this ob- 
servation and the almost complete bond-order parameter 
mismatch between the 4D crystal phases and the fluid it 
is reasonable to describe 4D as completely geometrically 
frustrated and 3D as only partially so. 

What about higher dimensions? Because the local 
bond-order mismatch is already very pronounced in 4D 
(see Sect. IIIBI) . it is hard to imagine that higher dimen- 
sions could show any more geometrical frustration than 
4D does. Once the fluid and crystal order distributions 
stop overlapping significantly, the fluid simply does not 
spontaneously form structures that easily anchor to the 
crystal interface. We therefore expect the coexistence 
fluid-crystal interfacial free energy to scale linearly with 

7™w> but onl y for d ^ 4 - 

Let us consider for a moment the impact of these obser- 
vations on 3D polydisperse systems. Increasing polydis- 
persity normally decreases Ava(p), as discussed above. 
The increase of 7f ° x x with polydispersity thus occurs 
because p cocx increases faster [83| than Ava(p° ock ) de- 
creases. The higher p coex and then lead to a higher 
free energy barrier and to a more rapid increase in non- 
equilibrium 7f_ x with supersaturation [84j . The curva- 
ture of the crystal nucleus appears to be a marginal con- 
tribution Hi]. 

E. Classical Nucleation Theory 

In the end, what do these results imply for the crys- 
tallization barrier? The geometrical interpretation for 
7f_ w above as well as the connection between 7f_ w and 
7f_ x suggested by geometrical frustration permit certain 
predictions. Because Ava(p) does not vanish in high 
dimensions and because p coox increases with dimension- 
ality, due to the increasing inefficiency of lattice pack- 
ings to fill space compared to simplex-based fluids, we 
expect 7f ° x x to grow with dimensionality. In the de- 
nominator of Eq. rT2] the crystal density increases with 
dimension, but the cell model provides a lower bound 
for the pressure contribution to 7f_ x that also scales lin- 
early with density and with a larger prefactor. Because 
the geometrical prefactor scales like d d l 2 , the nucleation 
barrier of monodisperse hard spheres thus increases with 
dimension. Higher-dimensional crystallization becomes 
ever rarer, which explains the surprising stability of su- 
persaturated monodisperse hard sphere fluids observed 
in simulation 0, [2l[ . 

III. CONCLUSION 

The modern understanding of geometrical frustration 
in hard sphere fluids considers how much space would 
need to be curved to allow for maximally dense simplex- 
based structures to form lattices 0, 0, S| ■ Though ulti- 
mately this lack of curvature is the reason why lattices 
in Euclidean space cannot be simple assemblies of sim- 



plexes, it lacks a direct dynamical mechanism to prevent 
crystallization. This study has made more precise the 
role of geometrical frustration in increasing the interfacial 
free energy between the fluid and the crystal in monodis- 
perse and polydisperse systems. We have argued that 
fluid-hard-wall interfacial free energy is an upper bound 
for the fluid-crystal interfacial equivalent, whose scaling 
behavior is easier to model. We have also argued that 
based on the poor overlap between the local fluid and 
crystal orders 7f_ x and 7f_ w should remain of the same 
order of magnitude in high dimension. Put together, 
these elements allow us to predict that in high dimensions 
the free energy barrier to crystallization is much larger 
than in 3D, and therefore crystallization is much rarer 
than in 3D. The crystal thus only marginally impacts a 
supersaturated fluid dynamics in high dimensions, unless 
a different type of phase transition arises [8(|. For the 
regime where that is not the case, higher dimensional 
spheres are an interesting model in which to study phe- 
nomena that are ambiguous in 3D, such as jamming, as 
we saw above, and glass formation [87j | . 

Why is crystallization then so common in 3D hard 
spheres? First, the crystal is sufficiently efficient at fill- 
ing space for the coexistence pressure to be relatively 
low. Second, the overlap of the fluid and crystal order- 
parameter distributions is significant in 3D, which re- 
sults in only a moderate geometrical frustration. Three- 
dimensional space is so small that all possible cluster or- 
ganizations, including the cubic crystal unit cells, are fre- 
quently observed in the fluid, which limits the extent of 
geometrical frustration. Bernal and many after him have 
indeed observ ed j ust about any small polyhedron in hard 
sphere fluids [88|]. Though the simplex-based fluid order 
is preferred, other ordering types are not far off and can 
easily be accommodated. Geometrical frustration in 3D 
monodisperse hard spheres is thus only partial, because 
even the crystal unit cell is sufficiently "liquid-like" . 

Finally, our study has shown that many microscopic 
details of the interfacial free energy of even the simplest 
of systems are mis- or incompletely understood. Any 
hope to get a grasp of and control homogeneous and 
heterogeneous crystallization in more complex systems is 
contingent upon having a better understanding of these 
fundamental issues. 
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APPENDIX A: LATTICE GENERATING 
MATRICES 



The canonical lattice representations are often given 
in a form that simplifies the notation [2^, Chap. 4], but 
does not necessarily allow for a convenient physical con- 
struction for simulation purposes. Moreover, the choice 
of generating lattice vectors should keep the size of the 
simulation box commensurate with the crystal unit cell 
to a minimum. Dd packings are checkerboard lattices, 
which are algorithmically simple to generate. A4 and A$ 
are dense packings with generating matrices 
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respectively. Eq is a cut through the E$ generalization 
of the diamond lattice, for which we use the generating 
matrix 
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(A3) 

From these generating matrices we can construct a unit 
cell commensurable with our hyper-rectangular simula- 
tion box. The sides of the unit cell Ui are obtained by 
linear combinations of the matrix' row vectors such that 
only one non-zero element remains. The lattice sites lo- 
cated within the unit cell borders then correspond to 
the particle positions within the unit cell. Following 
this recipe, our A4 unit cell has relative side dimensions 
I = (VlO, 1, \/2, 1) with n u — 8 particles in the unit cell, 
A 5 yields I = (Vl0, 1, \/3, Vl5) with n u = 120, and 
E 6 yields I = (3, 3, 3, \/3, y/3, y/S) with n u = 24. 



APPENDIX B: THIRD-ORDER INVARIANT 
POLYNOMIALS 



For reference, we provide two common 3D third-order 
rotationally-invariant polynomials w^(a, b, c) for 1=4 and 
1=6, expressed in inner products as outlined in Sec. II Bl 



The simple 4D I = 4 polynomial is also given. 
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APPENDIX C: STATISTICAL MEAN-FIELD 
THEORY OF JAMMING 

Using the notation, approach, and approximations of 
Ref. 221, we obtain the theoretical predictions of ?7mrj 
for higher dimensions. For mechanically stable config- 
urations, the bulk and contact contributions to the cu- 
mulative probability that the coordinates of all spheres j 
obey Tj j cos 9j > c is 



F> (c) w exp 



/ -2 d V*(c) 



-2dS*(c) 



\V d (l/r) stat 



MRJ 



1) S4-1/2 + Sd-2S oc 



with (see Table HV]) 

V*(c) = / e(c-r/cos6»)(ir (CI) 
S*(c) = J 6{r-l)S{c-r/ cos 6)dr (C2) 
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TABLE IV: Volume functions for the statistical mean-field jamming analysis of Ref. [22j] in higher dimensions. 
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TABLE V: Surface virial coefficients for hard spheres in low 
dimensions. 



A self-consistent solution to 



1/<RJ " 1 = " 



dc 



(C4) 



obtained numerically gives ?7mrj. 



The last result is obtained by noticing that the mirror 
image across the z — plane of a configuration is also a 
valid configuration and completes the lens formed by the 
addition of two spherical caps. The prefactor accounts 
for the double counting. 

To compute V£ (h), we use the spherical coordinates 

system Xd = rcos9 and \J x\ + ■ ■ ■ + x 2 d _ 1 = rsmO. 

Then the cutting plane x d — 1 — h intersects the sphere 
of radius 1 at angle 

Oh = arccos(l — h). 

Using the notation Sd-2 introduced in Section ITDl for the 
surface of a d— 2-dimensional sphere in d— 1-dimensional 
space, we compute 



Vpv(h) = S d -2 



r^ 1 sin^ 2 OdrdO 



APPENDIX D: HARD-SPHERE SURFACE 
VIRIAL COEFFICIENTS 



v ; cos d (6>) 



Using the notation of Ref. 63[ , we calculate the virial 
corrections to the interfacial free energy between a hard 
sphere fluid and a hard wall. The first two coefficients of 
the expansion 

Bm = 4/2 (Dl) 

S ra = 4/2-(24+4)/6 (D2) 

are given in Table fVl for d < 10. 

A convenient way to compute these coefficients uses 
the volume V^ ap (h) of a spherical cap of height h on the 
unit sphere x\ + ■ ■ ■ + x 2 d < 1. This cap is obtained as the 
portion lying above the plane Xd = 1 — h. 

We can then compute the quantities of Equations (|Dip 
and (|D2[) by use of the integrals 
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Recall that for J d {x) = J sm d (x)dx, we have 



J2m(x) 
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Using this information, we find 
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APPENDIX E: HEIGHT OF THE CAP CREATED 
BY THE INSERTION OF A WALL 

In this appendix, we compute the height h of the spher- 
ical cap created by inserting a wall in a densely packed 
simplex, which we use to approximate Ava- To simplify 
the computation, it is convenient to think of M. d as sit- 
ting in E d+1 as the hyperplane x\ + ■ ■ ■ + Xd+i = 
Let ei, . . . , &d+\ stand for the usual coordinate basis vec- 
tor. We put hard spheres centered at These 
spheres are all at a distance 1 from each other and form 
a (i-dimensional simplex. When inserting a wall tangent 
to the spheres centered at a/y/2, . . . , ed/V%, we cut a 
spherical cap from the sphere centered at e^+i/v^- We 
wish to compute the height of this cap. 

Let 



P d = 



ei 



dV2 



be the center of mass of the first d balls. Since Pd is at a 
distance yj 4^ from ed+i/V2, the wall being inserted is 

at a distance \J — \ from ed+i/V%- So the cap has 

height i-( A /^-i)=l 



Dividing this 

result by two for the two interfaces that are created by 
inserting a wall and taking the high-dimensional limit, 
we obtain 



l 

V2 



1 



Av A w -(2 - y/2) ~ 0.146447. 
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